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Abstract 

Over the last decade experimental studies have shown a large B isotope fractionation between materials carrying 
boron incorporated in trigonally and tetrahedrally coordinated sites, but the mechanisms responsible for producing the 
observed isotopic signatures are poorly known. In order to understand the boron isotope fractionation processes and 
to obtain a better interpretation of the experimental data and isotopic signatures observed in natural samples, we use 
first principles calculations based on density functional theory in conjunction with ab initio molecular dynamics and 
a new pseudofrequency analysis method to investigate the B isotope fractionation between B-bearing minerals (such 
as tourmaline and micas) and aqueous fluids containing H3BO3 and H4BO4 species. We confirm the experimental 
finding that the isotope fractionation is mainly driven by the coordination of the fractionating boron atoms and have 
found in addition that the strength of the produced isotopic signature is strongly correlated with the B-O bond length. 
We also demonstrate the ability of our computational scheme to predict the isotopic signatures of fluids at extreme 
pressures by showing the consistency of computed pressure-dependent /3 factors with the measured pressure shifts 
of the B-O vibrational frequencies of H3BO3 and H4BO4 in aqueous fluid. The comparison of the predicted with 
measured fractionation factors between boromuscovite and neutral fluid confirms the existence of the admixture of 
tetrahedral boron species in neutral fluid at high P and T found experimentally, which also explains the inconsistency 
between the various measurements on the tourmaline-mica system reported in the literature. Our investigation shows 
that the calculated equilibrium isotope fractionation factors have an accuracy comparable to the experiments and give 
unique and valuable insight into the processes governing the isotope fractionation mechanisms on the atomic scale. 
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1. Introduction 

The use of isotopes as geochemical tracers depends 
upon the existence of reliable P - T-dependent equi- 
librium isotope fractionation data between solids, fluids 
and melts. The common method used for determina- 
tion of such data is to conduct experiments, where the 
phases of interest are equilibrated at a range of different 
conditions and individually measured for their equilib- 
rium isotopic compositions. Recently, with the develop- 
ment of computational methods, software and increase 
in hardware performance it is also possible to simu- 
late and compute the isotope fractionation with ab initio 
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methods. These computational considerations comple- 
ment the experimental effort and provide information 
on the mechanisms governing the equilibrium isotope 
fractionation processes on the atomic scale. By estab- 
lishing an efficient computational approach for materi- 
als at high P and T and testing its reliability by com- 
puting Li isoto pe fractionation between minerals and 
aqueous fluids, Kowalski & JahnI (|201 ih have shown 
that ab initio methods can provide a reliable estima- 
tion of equilibrium isotope fractionation factors at an ac- 
curacy level comparable to experiments. Motivated by 
the encouraging results for the fractionation of lithium 
isotopes we applied our new method to study boron 
isotope fractionation. The main goals of the present 
study are to make theoretical predictions and obtain 
a better understanding of the B isotope fractionation 
process between tourmaline, B-bearing mica and flu- 
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ids at various pressures and temperatures, to compare 
these da ta to results f r om re c ent in situ experim ental 



studies (IWunder et all l2005t iMever et all l2008h and 



measurements of isotopic signatures in riatural samples 



jKlemme et al.L 1201 U iMarschalll 120051: iHervig et al. 



2002h and to investigate the underlying mechanisms 
driving the boron isotope fractionation processes be- 
tween the considered materials. 

With two stable isotopes '"B and "B, of relatively 
large mass difference of about 10%, boron isotopes 
strongly fractionate during geological processes, thus 
leading to natural ($'' B-variations ranging from -30 to 
+60 %o dBarthl Il993h . Therefore, B isotopes are ideal 
for the distinction of different geological environments 
and for quantifying mass transfer processes, e.g. in 
the range of subduction zones. First-order criteria driv- 
ing isotope fractionation in Earth materials are differ- 
ences in coordination and in the bonding environments 
of coexisting phases. The lighter isotope usually pref- 
erentially occupies the higher coordinated site, which is 
generally accompanied with a lon ger cation-anion bond 



WunderetaUl2011h . 



ige 

length and weaker b ond strength dSchauble et all 12009 



Tourmaline has an extensive chemical variability and 
is the most widespread borosilicate in various rocks over 
a large range of bulk compositions. It has a large P -T- 
stability which ranges from surface conditions to high 
pressures and temperatures of at least 7.0 GPa and about 
1000°C as determined experimentally for dravitic tour- 
maline ( Krosse, 1995 )■ Most tourmaline minerals con- 
tain 3 boron atoms per formula unit (pfu) with B in 
trigonal planar coordination (Bf-^'). High pressure com- 
bined with an Al-rich environment can lead to the for- 
mation of olenitic tourmaline with significant amounts 
of excess B substituting for Si at the tetrahedral site 
(gW-) 'pjjg highest amounts of up to 1.2 B''*^ pfu have 
been found in ole nitic tourmaline from Koralpe, Austria 
(lErtl et al.iri997l) . Al-rich tourmaline with up to 2.2 B'^^^ 
pfu was synthesized e xperimentally at 2.5 GPa, 600°C 
(ISchrever et al. ■ l2000h . The maximum possible amount 
of B^"^' in olenitic tourmaline is limited to three B'^^ pfu, 
due to structural and crystal chemical reasons. 

As tourmaline is not stable at basic pH 



([Morgan & Londonl Il989h . the dominant B-species of 
fluids coexisting with tourmaline is B(OH)3. Therefore, 
due to the absence of change in boron coordination, 
B -isotopic fractionation between B''*^-free tourmaline 
and fluid (A"B(tour-fl)) should be small. However, at 
0.2 GPa, the experimentally determined A"B(tour-fl)- 
values of -2.5 + 0.4 7m at 400°C and -0.4 + 0.4 %o at 
700°C (iMever et al .. 2008) suggest small but significant 
differences in the B-O(H) bond strength between 



tourmaline and the neutral fluids. Incorporation of B'^'*' 
into tourmaline is expected to significantly increase the 
fractionation of boron isotopes between the mineral and 
aqueous fluid (i.e. increase of |A"B|(tour-fl)). In this 
contribution we present calculated P - T-dependent 
data on B-isotope fractionation between B''*^-bearing 
tourmaline and fluid, which so far are not available 
from experiments. 

Despite of its low b oron content of up to maxi- 



mum values of 270 ppm (iDomanik et al.L Il993h . potas- 



sic white mica is probably the main host of boron 
in metasedimentary and metabasaltic blueschists and 
eclogites, because of its high modal abundance and 
B-incompatib ility in all oth er stable minerals of these 



rocks (.Brenan et al.L Il998h . During subduction the 



modal amount of mica continuously decreases by dehy- 
dration reactions and the chemistry of residual micas is 
shifted towards phengitic compositions. Phengitic mica 
has an extended stability, exte nding as deep a s 300 km 
within cold subduction zones (Schmidt, ' 1996*). In con- 
trast to most tourmalines, boron is tetrahedrally coordi- 
nated in mica, where it substitutes for aluminum. Due 
to the coordination change from mostly t hree-fold co 



ordinated in near-neutral fluids (Schmidt et al.. 2005) 



to B'"*^ in mica, B-isotopic fractionation between B- 
bearing mica and such fluids is much larger than for 
tourmaline - fluid. A"B(mica-fl)-values determined 



experimentally at P 
at 500°C and -7.1 d 



= 3.0 GPa, are -10.9 ± 1.3% 
0.5 %o at 700°C ('Wunder et al 



2005h . Such a strong B-isotope fractionation and its 



pronounced T-dependence, in combination with the 
continuous dehydration of micas during ongoing sub- 
duction and boron transport via fluids into mantle wedge 
regions of arc magma-formation, probably determines 
the boron variations and B-iso topic signatures in vol- 
canic arcs (IWunder et al.[l2005h . 

Using in situ Raman spect roscopic measuremeri ts of 
near-neutral B-bearing fluids, ISchmidt et al.l (l2005h ob- 
served a significant amount of 4-fold B-species at high 
P and T. The abundance of these species increases 
with temperature and pressure and at T - 800 K and 
P - 1.9 GPa there should be a considerable amount 
(15 - 30%) of Bl^l species in the fluid. Such a sig- 
nificant amount of tetrahedrally coordinated B-species 
in high temperature and pressure fluids should affect 
solid-fluid B-isotopic fractionation, which we investi- 
gate in our calculations. In the light of this, we also dis- 
cuss recently determined B-isotope data fi'om coexisting 
natur a l tourmalin e and B - bearing mica dKlemme et al.i 



20111; IMarschalll 120051; iHervig et al.L |2002'), which 
show slight inconsis tenc y with in situ measurements of 
Wunder et al.1 (l2005l) and lMever et"aD (l2008h . 
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Reliable ab initio computational methods to pre- 
dict isotope fractionation factors have been estab- 
lished recently. Several groups have proved that 
such calculations can contribute towards understand- 
ing geochemical mechani sms responsible for produc- 
tion of isotope signatures tPriesneii 1997t Yamaii et al 



2008 



| 200ll: ISchaubi3.l2004l; IPomagal-Goldrnan et al 
Hill & Schauble', '2008; 'Meheut et al.' '2009", '2007 
Schauble et al.- .2009: .Zeeba. .2009: .Hill et al.. 2010 



Rustad et all l2010al: iRustad et all l2010bl: IZeebd. 120091 
2OI0I : iKowalski & Jahn . 201 1 ). The majority of these 
works, however, concentrate on the computation of the 
stable isotope fractionation between various, mostly 
simple crystalline minerals, and the aqueous solu- 
tions are usually computed using an isolated clus- 
ter containin g fract i onatin g speci es and a h y dratio n 
shell (e.g. 'Zeebe' (2010', '2009'); 'Hill et al.' (2010'); 
Rustad et all (2()10a); .Domagal-Goldman et al (2008); 
£11 & Schaubld (l2008h~TRustad et all (l2010bh : IZeebe 
2005|)). However, aqueous solutions at high pres- 
sure and temperatur e must be computed w ith cau- 
tion as discussed in iKowalski & JahnI (1201 11) . This 
is because the distribution of cation coordination and 
cation-oxyge n bond lengths that affect the isotope 
fractionation jBigeleisen & Mavei . 1947 ) is driven by 
the dynar nics of the syst e m an d change unde r com 
pression dJahn & Wunder , I2OO9I: IWunderetall 1201 1 



Kowalski & Jahnl 201 ll). The only r ecent ab initio 
work, besides 'Kowal ski & JahnI ( 201 1 ). that accounts 



for the dynam ical effects on the i sotope fractionation 
in fluid is by Rustad & Bvlaskal (2007) who consid- 
ered boron equilibrium isotope fractionation between 
B(OH)3 and B(OH)4 species in aqueous solution. They 
performed ab initio molecular dynamics simulations of 
these fluids and attempted to use the vibrational density 
of states, derived through the Fourier transform of the 
velocity auto-correlation function, as an input for the 
calculation of the "B/'^B isotope fractionation coeffi- 
cient. However, the resulting fractionation factor a - 
0.86 happened to be much lower than the experimen- 
tal value of ff = 1.028. Interestingly, the discrepancy 
between experiment and theory is cured by quenching 
the selected configurations along the molecular dynam- 
ics trajectory and computing the harmonic frequencies. 
The fractionation factor derived using these frequencies 
exactly reproduces the experimental value. 

In our approach both solids and fluids are treated as 
extended systems by application of periodic boundary 
conditions in all three spatial directions, which is crucial 
to model high pressure materials. Large enough super- 
cells are chosen to avoid significant interaction between 
atoms and their periodic images, as well as to reduce the 



number of k-points and q-vectors for sampling the Bril- 
louin zones and the phonon spectra of the crystals, re- 
spectively (for a liquid or fluid, both the Brillouin zone 
and phonons are not defined). In our investigation we 
use cells of at least 7 A width in each spatial dimension. 
A representative statistical sampling of the fluid struc- 
ture is obtained by per forming Car-Parrinello m olecular 
dynamics simulations ( Car & Parrinellol 1985 ). For the 
calculation of the isotope fractionation factors in fluids, 
several random snapshots from the simulation runs are 
chosen. The force constants acting on the fractionating 
element and the resulting fractionation factors are then 
obtained for each ionic configuration, and the relevant 
fractionation factor for the boron species in the fluid 
is computed as an ave rage over the who l e set of con- 
sidered geometries. In Kowalski & Jahn ( 201 1 1) it was 
shown that in line with the i Bigeleisen & Maverl ( 1947 ) 
approximation, considering the force constants acting 
on the fractionating atom only leads to a satisfactory es- 
timation of the Li isotope fractionation factors for high 
temperature fluids and minerals. Here, we will show 
that approximating the vibrational spectrum by the three 
pseudofrequencies derived from the force constants al- 
lows for further improvement of the accuracy of the pre- 
dicted isotope fractionation factors, especially at lower 
temperatures. 

In this contribution we present the theoretical pre- 
diction of B isotope fractionation factors between B 
bearing aqueous fluids and solids, specifically tour- 
maline and B-muscovite, for which experimental data 
and measure ments on natural sani p les are available for 



comparison dWunder et al.L 120051; iMever et al.L 12008 



Klemme et al.Ll201 lb . We will show that the application 
of ab initio methods to B-bearing crystalline soUds and 
fluids not only provides unique insight into the mech- 
anisms driving equilibrium B -isotope fractionation on 
the atomic scale, but helps in proper interpretation of 
the data. 



2. Computational approach 

2.1. Theoretical model 

2.1.1. The single atom approximation: 

Bieeleisen &■ Maver approach 



Mass-dependent equilibrium isotope fractionation is 
driven by the change in molecular and crystalline vi- 
bration frequencies resulting from the different masses 
of the isotopes. The fractionation between species and 
an ideal monoatomic gas is called the /3 factor. In 
the harmonic approximation it is given by the formula 
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Table 1: The computed and measured vibrational frequencies for H3BO3 and H4B0^ molecules reported for two dilferent boron isotopes E'Ve'". 
Only the frequencies affected by the isotopic substitution are reported. The units are cm"' 



method 


H3BO3 




H4BO4 








theoretical: 














BLYP (this work) 


656/680 


1366/1408 


793/807 


870/886 


1051/1063 


1132/1154 


BLYP ' 


642/667 


1390/1437 


807/821 


916/947 


1020/1025 


1156/1179 


HF/6-31G*2 


736/764 


1562/1615 










HF/aug-cc-pVDZ* ^ 


-/754 


-/1531 










BP86/aug-cc-pVDZ* ^ 


-/655 


-/1412 










B3LYP/aug-cc-pVDZ* ^ 


-/684 


-/1447 










MP2/aug-cc-pVDZ* ^ 


-/684 


-/1435 










CCSD(T)/aug-cc-pVDZ* ^ 


-/684 


-/1438 










experimental: 
solution 


632/666 


1412/1454 




937/975 






vapour ^ 


674/700 


1429/1477 










solution ^ 


639/668 


1428/1490 




947/- 






vapour ^ 


666/692 


1415/1472 










vapour ^ 


675/701 


1426/1478 










solution and vapour ^ 


639-675/- 


1421-1450/- 




935-958/- 







(Il99lh . lOii (l2000l) : lLiu & Tosselll tOO± lAndrews & Burkholdedd 19921) . lOgden & Yound (Il988h . IZeebd (l2005h 

and references hereafter. 



2001): 



jBigeleisen & Maveii Il947l: lUrevll 19471: IChacko et al 



"dof t 

nu. 
- 



exp 



(m; - U*) 



1 - exp(-M,) 



1 - exp(-M*) 



(1) 



where u - hvilksT, h is the Planck constant, v, is the vi- 
brational frequency of the /-th degree of freedom, ks is 
the Boltzmann constant, Ndof is the number of degrees 
of freedom, which for being the number of atoms in 
the considered system (molecule, mineral or fluid) is 
equal to 3A^ - 5 for a diatomic molecule, 3A^ - 6 for 
multiatomic molecules and 3A^ for crystals, and a star 
symbol marks the heavier isotope. The fractionation 
factor between two substances A and B, qa-b is com- 
puted as the ratio of the relevant (3 factors, which is well 
approximated by the differences in the yS factors: 

(XA-B ^PaIPb, ^a-b = 10001n/?A - 10001n/3B[%o]. 

(2) 

The calculation of the /? factor requires only knowl- 
edge of the vibrational properties of the considered sys- 
tem computed for the two different isotopes. However, 
computation of the whole vibrational spectra of com- 
plex, multiparticle minerals or fluids requires substan- 
tial computational resources and is currently limited to 
systems containi ng a few dozens of atom s or less. In 
our recent work dKowalski&JahnlbOlli) we proposed 



to use an efficient method for computing the high tem- 
perature isotope fractionation factors between complex 
materials such as fluids and crystalline solids, which re- 
quires the knowledge of the force constants acting upon 
the fractionating element o nly. The P factor (Eg. [T) can 



be then approximated by (IBigeleisen & Maverl 11947 
■Kowalski & Jahn. 2011): 



/3- 



^ ?4 



i=I 



1 + 



Am 



mm* 24klT^ . , 



where A, are the force constants acting on the isotopic 
atom in the three perpendicular spatial directions (x, y 
and z). Am - m* - m, where m and m* are the masses of 
the lighter and heavier isotopes of the fractionating ele- 
ment. As the computation of the /3 factors from formula 
|3] requires the knowledge of properties of the fraction- 
ating element only we will call such an approach the 
single atom approximation throughout the paper The 
validity criteria restricts the usage of the formula to fre- 
quencies v[cm~'l < 1.39 r[K 1 (assu ming u < 2, see 
Fig. 1 of Bigeleisen & Maver ( 1947h '). We are inter- 
ested in temperature range 800-1000 K. The highest vi- 
brational frequency of the modes involving movement 
of B atoms for H3BO3 is ~ 1400 cm " and of H4B04' 
is < 1200 cm ■ (Table [D. In case of H3BO3, the sin- 
gle atom approximation may produce an error of 2.2 %o 
in the /3 factor at T = 800 K. The relevant error for 
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Table 2: WQ()(/3— 1) factors for isolated H3BO3 and H4BOJ m olecules at T = 300 K obtained using three methods: (1) the full frequency spectrum 
and Equation[T] (2) the single atom approximation of Kowals ki & Jah j )201 ih . (3) the single atom approximation with the pseudofrequencies and 
equation^ The units are %o. The A/3//Jbm is the check of condition given by Eq. llSl and indicates the improvement of the method (3) over method 
(2) expressed in %. 



T(K) 


meth. (1) 


meth. (2) 


meth. (3) 


^IPbm [%] 


meth. (1) 


meth. (2) 


meth. (3) 


(^PIPbm [%] 




H3BO3 








H4BO4 








300 


211.3 


283.1 


225.1 


20.5 


167.3 


203.4 


175.2 


13.9 


600 


64.1 


70.8 


65.6 


7.3 


48.1 


50.9 


48.6 


4.5 


800 


37.6 


39.8 


38.1 


4.3 


27.8 


28.6 


27.9 


2.5 


1000 


24.6 


25.4 


24.7 


2.8 


18.0 


18.3 


18.0 


1.7 




Dravite 








Boromuscovite IM 








300 


202.8 


267.7 


216.0 


19.3 


139.1 


160.8 


142.2 


11.6 


600 


60.7 


66.9 


62.4 


6.7 


38.6 


40.2 


38.8 


3.5 


800 


35.5 


37.6 


36.1 


4.0 


22.1 


22.6 


22.2 


1.8 


1000 


23.2 


24.1 


23.5 


2.5 


14.3 


14.5 


14.3 


1.4 



H4BO4' is 0.8 %o (Table 121). A further improvement to 
the method is therefore desired. 

2.7.2. The single atom approximation with pseudofre- 
quencies: our improvement 
We will show that the error of the single atom ap- 
proximation can be substantially reduced if one uses the 
three frequencies v,- derived from the force constants 
acting on the fractionating element (v^ = AijAn^m). 
We call them "pseudofrequencies ", and compute the fi 
factors using formula [1] In the following we present 
the formal justification of such an approach. Accord- 
ing tolgigeleisen & May er ( 1 947), equation[T]for smal l 
Am,- - Ui - u* reduces to (IBigeleisen & Mavei (Il947h . 
Eq. 11a): 

/II 1 \ 

>8=i+yu--+^^^^«'- (4) 

j-^\2 Ui exp(M,)-l/ 

The Taylor expansion of the function appearing under 
the summation sign is: 

G(u) = i - i + ^ 



3 5 
u u u 



2 u exp(M,) - 1 

„7 



12 720 ■ 30240 ' 1209600 "" 
When we consider just the first term of the expansion 
the p factor is: 

N,i„f 



^=i+;^-^A«,. 



Ndof A 2 ^<iof 2 *2 

Am -t— ; m - M 



i= 1 /= 1 



24 ' 



(6) 



Table 3: The three pseudofrequencies of the B atom derived for two 
selected molecules and two crystalUne soUds. The units are cm"' . 



species 


H3BO3 


644/675 


1129/1184 


1130/1185 


dravite 


618/648 


1085/1138 


1120/1174 


H4BO4 


809/848 


845/886 


872/915 


boromuscovite IM 


722/757 


723/758 


801/840 



which is exactly equation |3] 

Let us consider the Taylor expansion of the different 
estimations of fi factors. Equation |4] then reads: 



= 1 + 



/ ,,3 
^ I 12 720 

/=1 ^ 



Am,-. (7) 



The tei geleisen & Mava ( 1947 ) approximation given 
by equation[3]reads: 



(8) 



and the proposed approximation based on pseudofre- 
quencies is: 



P ^ Ppseudo 



3 / - 
Y"! I M, 



=1 ^ 



12 720 



Aui. 



(9) 



In the above equation m, = hvi/ksT. Next we check how 
this ap proximation compares to the Bi geleisen & Mavei 
(119471) approximation given by equations [3] and [S] In 
order to make a comparison we derive the differences 
between the two approximate expressions /3bm, Ppseudo 
and the exact one fiexact (Eq. Q. In the case of the 
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1000 



600 
T(K) 



1000 



Figure 1: fi factors for isolated H3BO3 and H4B0^ molecules as well as dravite and boromuscovite IM crystallin e solids. The lines represen t 
the results obtained using: (1) the full frequency spectrum and Eq. [T](solid), (2) the single atom approximation of lBigeleisen 
iKow alski & Jahn ( 2011J, Eq.|3](dotted) and (3) the method described in the text with the pseudofrequencies and Eq. ^(dashed lines). 



Bigeleisen & Maveii (Il947h approximation we have: 



Zili _ V ill _ _i_ 
12 ZjIi2 720 



12 

V I 



Am,- 



Am,- 



11^720 

and having from equations |3]and[8]that 

Njof 3 _ 

— Am,- — > — Am,- 
12 ^12 



(10) 



(11) 



in the case of the proposed approximation we get; 

3 / .-,3 N 



^Ppseudo — Ppseudo Pexact — ^ ] |^ 720 



Am/ 



^720 



1^720 



^\ 720 



Am,- = 



Am,- + AySfiM 



(12) 



Because relation G (u) < -f^ holds for any m (see 

Bigeleisen & Maveii(ll947l) . Fig. 1), the function 



720 



Am,- = Yj - ^) < (13) 



and ls.lipseudo < ^Pbm- On the other hand the expres- 
sion for A/3pseudo is given by a difference of the higher 
order terms of the Taylor expansions of the two expres- 
sions for the /3 factor. In the considered cases the values 
of pseudofrequencies are similar to the real frequencies 
that are affected upon B isotope substitution in a given 
B-bearing system. This can be seen by comparing the 
pseudofrequencies computed for the selected cases of 
B-bearing molecules and crystalline solids considered 
here and reported in Table [3] with the real frequencies 
given in Table [1] This indicates that the two terms of 
opposite signs in Eq. [T2] should be similar in value and 
cancel out to a great extent, so \^i^pseudo\ « ^Pbm- 
Therefore the approach proposed here to compute the 
P factor based on pseudofrequencies and Eq. [T] should 
give a better approximation to the exact p factors than 
equation[3] wh ich we will show in sectionl3.1l The dif- 
ference to the iBigeleisen & Maver ( 1947h approxima- 



Table 4: The lattice parameters of the investigated B-bearing crystalline solids. Nawms is the number of atoms in the modeled supercell. 



a (A) 
b(A) 
c(A) 

an 

pn 

^ atoms 



boromuscovite IM boromuscovite 2M1 dravite olenite 



10.204 

8.788 

10.076 

90 
101.23 

90 

84 



10.180 
8.822 
19.8189 

90 
95.62 

90 

168 



15.945 
15.945 
7.210 

90 

90 

90 

163 



15.5996 
15.5996 
7.0224 

90 

90 

120 

162 



References: ' iLiang et alj (il995i) . ^Ertl et al.i (i2010i) . ^Marler et alJ (12002) 



tion is given by: 



(14) 



1=1 



and can be easily computed for any considered system. 
We assume that the pseudofrequency -based approach to 
the computation of fi factors is applicable if it is just a 
correction to equation[3] i.e. when: 



12 

/=1 



(15) 



We also note that the prop osed approach s atisfies the 
Redlich-Teller product rule dRedlichl 1 1935b when the 
mass of considered isotope m is much smaller than the 
mass of the whole considered system M, namely, 



U: I M m ^ 



nil. 



M m* 



(4r=n-- (16) 



We notice that we could force the strict conservation of 
the Redlich-Teller product by just adjusting the ratios of 
Uilu*. However, such a modification would not preserve 
relation [TT] and the pseudofrequency approach would 
not recover exactly the high temperature limit (Eq. |3]l 
of the exact solution (Eq. [U, which is a more important 
constraint to fulfill strictly by the proposed approxima- 
tion. 

2.2. Representation of solids 

In this paper we investigate the boron isotope frac- 
tionation between dravite, olenite and boromuscovite 
minerals and aqueous fluids. The solids were rep- 
resented by large cells containing at least 84 atoms. 
The number of atoms used in the crystal calculations 
together with the lattice parameters of modeled crys- 
tals are summarized in Table |4] The lattice parame- 
ters and chemical compositions of the modeled crys- 
talline solids are the experimental values measured at 



ambient conditions found in the literature. Dravite 
is the crystalline solid which was used in the experi- 
ments on tourmaline bv iMever et al.l (120081) . The chem- 
ical composition of the supercell used in the inves- 
tigation is Na3M g9Ali8(Sii8054)(B'^ ^'03)9(OH)i2 with 
structural data of Marler et al.l (120021) . Olenite can con- 
tain B in both trigonal a nd tetragonal sit es. The mod- 



eled structure is diat of lErtlet all (I2OIOI) . The chem 



ical composition of the unit cell used in the investi- 
gation is NaAl3Al6(Si4B['*^Oi8)(B[3]03)3(OH)30. For 
boromuscovite, th e IM and 2M1 crystal structures of 
Liang et al.l (Il995h were used. In the i sotope fractiona- 
tion experiments of lWunder et al.l (12005.) boromuscovite 
forms two polytypes, IM and 2M1, with relative abun- 
dances of 10% and 90% respectively. In boromus- 
covite B occupies the 4-fold coordinated site occupied 
mainly by Si atoms. The constructed model constitutes 
a 2x1x1 supercell of elementary chemical composition 
KAl2(BmSi30io)(OH)2. 



2.3. Representation of aqueous solution 

The aqueous solution was represented by a periodi- 
cally repeated box containing up to 64 water molecules 
and one H3BO3 or H4BO4 molecule. The pressure 
and temperature conditions w ere chosen to be clo se to 
the exper imental conditions of lWunder et alJ (120051) and 
Mever et al.. (|2P08). The pressure of the aqueous solu- 
tion for a given temperature and volu me was calculated 
according to the equation of state of IWagner & Pruss 
(.2002 ). The ab initio molecular dynamics simula- 
tions (AIMD) of aqueous fluids were performed for 
fixed temperature and volume using the Car-Parrinello 
scheme (Car & Parrinellp, 1985). The temperature dur- 
ing each ru n was controlled by a Nose-Hoover chain 
thermostat ('Nose & Kleliii 1 19831: iHooveii Il985b . For 
each T - V conditions at least lOps long trajectories 
were generated with an integration step of 0.12 fs. 
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Figure 2: Left panel: yS factors for H3BO3 (solid line) and H4B0^ (dashed lines) in the gas phase. Right panel: fractionation factors b etween 
H3BO 3 and H4BO^ in gas phase. The thick lines represent our results, while the thin lines represent the results obtained using frequencies of lZeebd 

0oo3) . 



2.4. Computational technique 

The calculations of pseudofrequencies and fi fac- 
tors for solids and aqueous solutions were performed 
by applying density functional theory (DFT) methods, 
which are currently the most efficient methods allowing 
for treating extended many particle systems quantum- 
mechanically. We used the planewave DFT code CPMD 
dMarx & Hutteii l2000l) . which is especially suited for 
ab initio simulations of fluids, the BLYP exchange- 



correlation functional jBeckel Il988t iLee et al.L Il988l) 



and norm-conserving Goedecker pse udopotentials for 



the d escription of the core electrons jGoedecker et al. 
1996h . One advantage of using the BLYP func 



tional is that it usually gives harmonic frequencies 
that most closely resemble the o bserved frequencies 
of benchmark ch emical systems (' pinlev & Stephens , 



I995I: lAlecu et all I2OIO.) Q. The energy cut-off for the 



plane wave basis set was 70Ryd for geometry relax- 
ations and molecular dynamics simulations and 140 Ryd 
for computation of vibrational frequencies. Periodic 
boundary conditions were applied for both crystalline 
solids and aqueous solutions to preserve the continuity 
of the media. 



' If the computed harmonic frequencies are closer to the observed, 
anharmonic frequencies than to the real harmonic frequencies of the 
computed system the resulted fractionation factor fi (Eq. Q computed 
from these frequencies should be closer to the real fractionation factor 
than the j6 computed on a set of accurate harmonic frequencies. 



The force constants and frequencies needed for the 
computation of the ji factors were computed using the 
finite displacement scheme. Before performing the cal- 
culations of the crystal structures all atomic positions 
were relaxed to the equilibrium positions to minimize 
the forces acting on the atoms. We note that to com- 
pute the fi factors for crystals one formally should ac- 
count for phonon dispersion. Here we use large su- 

percells and restrict our calculations to a single phonon 

\ 1 \ — — I 

wave-vector (F). Schauble (2011) has shown recently 

for ^^Mg/^'*Mg fractionation in Mg -bearing minerals 
that supercells containing more than 20 atoms are suffi- 
cient to get very accurate fi factors even at T = 300 K 
(eiTor of 0.1 %c). At r = 1000 K the error is in the 
order of 0.01 %o. The accuracy of the high tempera- 
ture isotope fractionation factors computed on a sin- 
gle phonon wave-ve ctor is also demonstrate d for iron- 
bearing minerals by iBlanchard et al. I (l2009h and con- 
firmed with good agreement of the predicted with the 
measured Li isotope fractionation factors between stau- 
rolite, spodumene, micas and aqueous fluid presented in 
our previous work (iKowalski & JahnLl201 ih . 



Prior to the computation of the force constants and 
frequencies of boron atoms in the fluids the positions of 
all the atoms constituting the boron-carrying molecule 
(H3BO3 or H4BO4) were relaxed to the equilibrium po- 
sitions, while all other atomic positions remained un- 
changed. The full normal mode analyzes were per- 
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Figure 3: Left panel: /J factors for H3BO3 (solid line) and H4B0^ (dashed lines) in aqueous solution. Right panel: fr actionation factors between 
H3BO3 and H4BO^ in aqueous solution. The thick lines represent our results, while the thin lines rep resent the results oljSanchez-Valle et al.l feOOSi) 
obtained using harmonic freque ncies (their Table 2). The dotted line represents the corrected Sanch ez- Valle et ajl j200a) results. The correction is 
made by comparing the work of lRustad et all J2010bl) and the correction derivation procedure is discussed in the text. 



Table 5: The 1000(/3- 1), ff = /3h,B03 /jSh4B0j and A = 1000(ln/?H3BO3 -In/JHiBO-) factors computed for isolated H3BO3 andH4BO;, atT = 300 K 
bv lZeebdtlOOl using different basis sets and our result obtained using the full normal mode spectrum and equation[T] The units are %o. 





H3BO3 


H4BO4 


« - /5h,B03//5h4BO- 


A = 1000(ln;6H3BO3 - lnj8H4BO-) 


6-31+G(d) 


216.6 


174.4 


1.0359 


35.3 


6-311+G(d,p) 


215.0 


170.5 


1.0380 


37.3 


our work 


211.3 


167.3 


1.0377 


37.0 



formed using the same method, but displacing all the 
atoms constituting the considered system. In the latter 
case the frequencies were obtained th rough the diago- 
nalization of the full dynamical matrix (ISchaubld.120041) 
as implemented in CPMD code. The effect of the vari- 
ous approximations on the derived fractionation factors 
was studied by additional computations of H3BO3 and 
H4BO4 isolated clusters. For that purpose we used a 
large, isolated simulation box with a cell length of 16 A, 
forcing the charge density to be zero at the boundary, 
as implemented in CPMD code. In order to compute 
the /3 factors of boron species in the aqueous fluid we 
apply the same method as in ou r recent work on Li 
isotopes dKowalski & Jahnl |201 ih . with the exception 
that we use the pseudofrequencies, i.e. the frequen- 
cies obtained from the three force constants acting on 
the fractionating element, and formula[T]for calculation 
of j3 factors, as discussed in section 12.1.21 In order to 
fully account for the spatial continuity of the fluid and 



its dynamical motion we produced 10 ps long molecular 
dynamics trajectories of systems consisting of 64 H2O 
molecules and one H3BO3 or H4BO^ molecule for dif- 
ferent T = lOOOK, 800 K and 600 K and pressure of 
0.5 GPa, \ yhich closely resemb les t he experimental con - 
ditions of IWunder eTal] (l2005h and lMever eTal] (l2008h . 
The corresponding simulation box length is 13.75 A at 
T - 1000 K. The /3 factors were computed on the ionic 
configuration snapshots extracted uniformly in 0.1 ps 
intervals along the molecular dynamics trajectories. 



2.5. Error estimation technique 

The errors in the computed value of the (/? - 1) and 
A fractionation factors were estimated from an aver- 
age error of vibrational freque ncies computed usin g 
the chosen DFT method . Finlev & StephensI (1995 ). 
Menconi & Tozej (|2002|) and lAlecu et alj (I2OI0I) esti- 
mated the errors made in calculations of vibrational fre- 
quencies of small molecules using different DFT func- 
tionals. According to these works the BLYP functional 
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systematically overestimates the harmonic frequencies 
by ~ 3.5 %, with a deviation from the mean offset of 
~ 1 %. Therefore, we expect that using BLYP func- 
tional the (J3 - 1) and A values are systematically over- 
estimated by 7 % and that in addition there is a 2 % 
error in derived (J3 - I) factors. Similar errors result 
from using other functionals or even more sophisticated 
and tim e consuming post-Hartree-Fock meth ods such as 



MP2 (.Finlev & Stephenslll995HAlecu et al.Ll201Q) 



3. Results and discussion 

3.1. Test of the computational method 

First, we illustrate the performance of the approxima- 
tion proposed in section |2. 1.21 bv computing the y6 fac- 
tors for the isolated H3BO3 and H4BO4 molecules and 
selected crystalline solids. In Figure[T|we present three 
sets of calculations of j3 factors: (1) the "exact" result 
obtained from a full normal mode an alysis and formula 
[Tj (2) the results obtained applying iKowalsk i & Jahn 
( 1201 lb method based on Eq. [3] (3) the results obtained 
using pseudofrequencies computed for the fractionating 
element and Eq. [T]for the estimation of the (3 factor The 
numerical values for selected temperatures are reported 
in Table |2l Approach (3) results in much better agree- 
ment with the "exact" result. For H3BO3, the fi factor is 
overestimated by only 0.5 %c and 1 .5%c for temperatures 
of 800 K and 600 K respectively. Applying method (2), 
the error is more pronounced, 2.2%o and 6.7%o respec- 
tively. In the case of molecular H4BO4, the errors using 
method (3) for the same temperatures are only 0.1 %o 
and 0.5 %c respectively. The same behavior is shown 
for dravite and boromuscovite crystalline solids that 
contain boron in the coordinative arrangement that re- 
semble the configurations of aforementioned B-bearing 
molecules. For T > 600 K the proposed method repre- 
sents only a few percent correction to the approximation 
given by equation (O, so the relation ( fTSl l is satisfied. It 
is evident that for B-bearing materials considered here 
the improvement made by using the pseudofrequencies 
based approach is substantial. It corrects for about 75% 
of error of the £igeleisen & Mayer ( 1947) approxima- 
tion (Eq. O. However, the question of general applica- 
bility of the proposed method to other isotopic systems 
would require careful testing on a large set of materials, 
which is well beyond the scope of the current paper. 

3.2. B isotope fractionation in gas and fluid phases 
3.2.1. B isotope fractionation between H3BO3 and 

H4BO4 in the gas phase 
In a first step of our investigation of boron-rich aque- 
ous fluids we derived the full frequency spectra of 



molecules in the gas phase (isolated molecules). The 
relevant /3 factors were computed using Equation [1] 
These studies were performed in o rder to compa re our 
results with the published values of lZeebd (l2005h . both 
computed using the same DFT BLYP functional. In Ta- 
ble [T] we report the computed frequencies that are af- 
fected by the different B isotope substitutions along with 
other theoretical estimations and experimental measure- 
ments. The computed frequencies are in good agree- 
ment with earlier theoretical predictions and show simi- 
lar agreement with the experimental measurements. The 
results in terms of computed /3 factors for the two con- 
sidered species are reported in Figure|2] where we com- 
pare our re sults w i th the values computed using fre- 
quencies of IZeebd (l2005h . The comparison of the two 
sets of calculations reveals that our yS factors for both 
species are smaller by ~ 1 %o at 600 K - 1000 K than 



the values of Zeebd (120051) . However, the difference be- 
tween j3 factors of H3BO3 and H4BO4 remains nearly 
identical in both sets of calculations and the agree- 
ment is nearly perfect for higher temperatures. We 
no te, that for th e comparison we used the frequencies 
of IZeebd (l2005h computed using 6-31-i-G(d) basis set, 
as only these are provided by the authors. /3 and a fac- 
tors obtained at T - 300 K using a more extended 6- 
311 -i-G(d,p) basis set indicate that the /3 factors using 
6-31-i-G(d) basis set are not fully converged. In par- 
ticular, the (J3 - I) factor of H4BO4 computed with 6- 
311-i-G(d,p) basis set is 3.9 %c smaller than the one de- 
rived using 6-31-HG(d). In Table |5] we compare these 
results with the results of our calculation. It is clearly 
seen that for lower temperatures such as T = 300 K 
the values computed with 6-311-i-G(d,p) basis set are in 
better agreement with our results indicating that plane- 
wave based DFT approach we use provides adequate vi- 
brational frequencies and resulting isotope fractionation 
factors. 

3.2.2. B isotope fractionation between H3BO3 and 
H4BO4 in aqueous fluid 

In order to obtain the temperature dependent /3 factor 
for aqueous fluids we fitted the function 1 +A/T^ + B/T'^ 
to the computed values using the least squares mini- 
mization procedure. The computed /3 values for H3BO3 
in fluid are: 1.02366 + 0.00012, 1.03624 + 0.00018 and 
1.06262+0.00010andforH4BO4 in fluid are: 1.01745 + 
0.00005, 1.02650 + 0.00010 and 1.04597 +0.00015, for 
the temperatures of lOOOK, 800 K and 600 K respec- 
tively. The resulting temperature dependent /3 factor for 
H3BO3 is j6 = 1 + 2.416 ■ lOVr^ - 5.823 ■ lOV^"* and 
for H4BO4 is /? = 1 + 1.772 ■ 10* /T^ - 4.234 ■ lO^T'^. 

The results for H3BO3 and H4BO4 in aqueous so- 
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Figure 4: Left panel: pressure dependence of the yS factor of neutral fluid (H3BO3) at T = 1000 K. Dashed line is the lineal' regression fit to 
the calculated values (points): iOQO(j3 - 1) = 23.6 + 0.28 P(GPa); Right panel: The computed change of the jS factor with pressure (symb ols) in 
compa rison to the increase in the fS factor derived from the frequency shifts of the 666 cm"' and 1454 cm"' lines measured bv lSanchez-Valle et alJ 
420051) (dashed hne). The coinparison is made assuining that (fi-l)~v^. 



lutions are shown in Figure [3] As was observed for 
the isolated molecules, the /3 factor of HaBOs-bearing 
fluid is substantially larger than the one for the H4BO4 . 
This can be understood in terms of the substantial dif- 
ference in the B-O bond lengths exhibited by the two 
considered species. In case of isolated molecules our 
calculations indicate a B-O bond length of 1.40 A for 
H3BO3 and 1.51 A for H4BO4. We compared our 13 fac 



tors w ith the values computed by ISanchez-Valle et alJ 12007; Liu & Tossell , 2005; Zeebe, 2005). Rustad 



also plotted in Figure [3] It is now very consistent with 
our prediction. 

3.2.3. Discussion of computational errors 

Most previous computational studies of boron iso- 
tope fractionation in aqueous solutions concentrate on 

the computation of the isotope fractionation at ambi- 

i 1 1 1 1 

ent c onditions (Rustad et al, 2010b; Rustad & Bvla skal 

etall 



(l2005h . which were derived by the combination of force 
field methods and experimental data to derive accurate 
vibrational frequencies. For H3BO3 we got a nearly 
identical result. In case of H4BO4 our calculation pre- 
dicts a valu e which is l ower by 2 - 4 %o. However, 
Rustad et all (l2010bl) and iRustad & Bvlaskal (|2007) re- 
vealed the improper assignment of a major fractionat- 
ing vibrational mod e of H4BO4 in the force field by 



Sanchez- Valle et al 



12005 ). This leads to the under- 
estimation of the fract ionation factor between aqu eous 
H3BO3 and H4BO; bv lSanchez-VaUe et al.1 (l2005h . As- 
suming that a oc T^^ and having the d i fference be- 
tween BLYP calculations of Rustad et a 1 (1201 Obi) and 
Sanchez- Valle et all (l2005h of Aa = 16.4 %o at T = 



300 K, tiie value reported bv lSanchez-VaUe et al.1 (120051) 
should be underestimated by Aa = 16.4 ■ 2(300/ r)^ %o, 
which results in Aa ~ 1.5 %o atT - 1000 K. Corrected 



in such as way result of ISanchez- Valle et al.l (l2005l) is 



(l2010bb performed detailed analysis of impact of the 
chosen computational method (HF, MP2, different DFT 
functionals) and size of the basis set on the calculated 
fractionation factors between H3BO3 andH4B04. They 
found that DFT methods are not performing well for 
the borate system and concluded that DFT "is of lim- 
ited usefulness in chemically accurate p redictions o, 
isotope fractionation in aqueous systems " (IRustad et a 
2010b,) . The empirically derived error of the derived 
fractionation factor is of the order of 5 - 10 %c for a 
total fractionation of ~ 30 %c. We note that this is ex- 
pected and clearly visible if we apply the error estima- 
tion procedure outlined in section 12.51 For instance, 
at room temperature the derived beta factors using the 
BLYP functional are 213.6 %o and 173.3 %o respectively 
(IRustad et all (l2010bl) . Table 2). This gives a fraction- 
ation factor of 1.0343. Following our error estimation 
scheme, the absolute error of the fractionation factor 
is 10.5 %c, and the properly reported computed value 
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Figure 5: Left panel: pressure dependence of the yS factor of strongly basic fluid (H4B0^) at T = 1000 K. Dashed line is the regression fit to 
the calculated values (points): 1000(/? — 1) = 17.15 + 0.754 P - 0.027 P-, where pressure is given in GPa; Right panel: The computed change 
of the factor with pressu re in comparison to the increase in the fS factor derived from the frequency shift of the 975 cm" ' fine measured by 
ISanchez-Valle et alJilOO^ (dashed line). The compaiison is made assuming that (jS - 1) ~ v . 



is a - 1.034 + 0.011. When one corrects for the sys- 
tematic error of 7% and assumes 2% of statistical error 
on yS factors, then the value of a decreases and the er- 
ror is slightly smaller, i.e. a - 1.032 + 0.008. This is 
in good agreem ent with the experimental data reported 



Rustad et a 1 (l2010b ) and explains the spread of the 



m 

values computed using different methods and reported 
in that paper. It is very difficult to get the fractionation 
factors for ambient conditions, as the fractionation fac- 
tor is often just a small fraction of the relevant (J3 - 1) 
factors, (a - 1) ~ 0.15(/3 - 1) in the considered case. 
Assuming that (a - I) = 0.15(yS - 1), a 2% error in 
the (J3 - 1) factors leads to an absolute error in (a - 1) 
of 0.04(/? - 1) = 0.04(0' - 1)/0.15 ~ 0.27(0' - 1), i.e. 
~ 27% of relative error in the derived fractionation fac- 
tor (o- - 1). On the other hand, we note that such a big 
error is not substantially larger than the uncertainties in 
the experimental data reported by iRustad et al (1201 Obi) 
in their Figure 2. Thus, the case of boron fractionation 
in aqueous fluid at ambient conditions does not neces- 
sarily show the limited usefulness of DFT in the predic- 
tion of isotope fractionation factors, but only reflects the 
fact that precise estimation or measurement of the B iso- 
topes fractionation factors at ambient condition requires 
unprecedented accuracy of both experimental or com- 
putational techniques. For instance, in order to get the 
value of (a - 1) with a relative error of 5% (at ambient 



conditions) one needs to estimate the (yS - 1) factors or 
measure relevant quantities with precision of less than 
1%. At higher temperatures the situation is different. 
Looking just at the fractionation between H3BO3 and 
H4BO4 in the gas phase or the aqueous solution one 
can see that for T > 600 K the fractionation factor be- 
tween the two substances, (a - 1), is at least 25% of the 
(yS - 1) factor. This results in smafler 0.04/0.25 ~ 16% 
for r = 600K and 0.04/0.36 ~ 11% for T = 1000 K 
relative error, which is acceptable in our calculations. 
Nevertheless, this case shows the importance of proper 
error estimation on the computed fractionation factors. 
Such an estimation is usually omitted or not provided 
explicitly, which can lead to wrong conclusions when 
the theoretical prediction is confronted with the mea- 
sured data. 

3.2.4. Pressure dependence of the fluid fractionation 
factor 



In our recent paper dKowalski & Jahnl 1201 Ih we 
have shown that due to compression the p factor of 
Li in aqueous fluid increases with increase in pres- 
sure (for P > 2 GPa). The same should happen for 
H3BO3 and H4BO4 aqueous fluid as the vibrational 
frequencies of boron speci es in aqueous fluid increase 



2005 



with increase in pre ssure jSanchez-Valle et al. 
Schnjidt et all 2005 ). Having the experimental data 



we checked whether the derived pressure-dependent /? 
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factors are consistent with the pressure shifts of vi- 
brational fregu£ncies_of_coimdered^ species mea- 
sured by |SM£hez::WkeLalJ (|2005|). For that purpose 
we performed a set of calculations using supercells 
containing 8 water molecules and the relevant boron 
species. We note that in line with our previous re- 
sults for Li (iKowalski & Jahnl 1201 ih . the obtained val- 
ues of (y8 - 1) at P = 0.5 GPa are within 0.1 %o in 
agreement with the values obtained for supercells con- 
taining 64 water molecules. The results are given in 
Figures |4] and |5] The computed (J3 - I) values for 
H3BO3 fluid show a linear dependence in pressure, 
(J3 - I) ^ 23.60 + 0.28 P(GPa) %o. This i s exp ected 
as (jS - 1) oc v2 ~ vl + 2voAv ("Schauble, '2004") and 
Av is a linear functio n of pre ssure dSanchez-Valle et al.L 
20051: ISchmidt et ail l2005h . In case of H4BO- the 
pressure-dependence is linear up to f ~ 2-3 GPa 
and it becomes less steep at higher pressures. In order 
to quantitatively check the consistency of our predic- 
tion with the mea sured vibrational frequency shifts of 
Sanchez- VaUe et al. (2005 ) we derived the relative shifts 



in the (/? - 1) factor assuming that (J3 - 1) oc and the 
measured pressure dependence of the frequency shifts: 
Av = 2.15 cm-' ■ P(GPa) and Av = 3.50 cm"' ■ P(GPa) 
for 1454 cm ' and 666 cm ' vibrational frequencies of 
H3BO3 and Av = 6.47 cm ' ■ P(GPa) for the 975 cm " 
vibrational frequency of H4BO4. The chosen vibra- 
tional frequencies are the ones affected by the diff'erent 
B isotope substitution. Our predicted shift of (J3 - 1) 
matches well the shifts derived from the measured fre- 
quency shifts. Such a good agreement with the ex- 
perimental data validates further our computational ap- 
proach and shows that ab initio calculations can be suc- 
cessfully used in derivation of the pressure dependence 
of the fractionation factors and pressure-induced vibra- 
tional frequencies shifts. Moreover, first principles cal- 
culations can be useful in extrapolation of the experi- 
mental values for j6 and Av to more extreme conditions, 
which otherwise are extremely difficult to reach by ex- 
perimental techniques. 

3.3. Fluid-mineral fractionation 

Next we present the results of the fractionation be- 
tween boron bearing fluids and minerals such as dravite, 
olenite and boromuscovite. The aim of these studies 
is to investigate the mechanisms driving the fractiona- 
tion process, the role of coordination and the B-O bond 
length. Below we discuss each case separately. 

3.3.1. Tourmaline-neutral fluid 
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Figure 6: The fractionation factors between tourmaline and aqueous 
fluid. The solid Une represents our prediction for the fractionation 
between dravite and H3BO3 neutral fluid and the shadowed region 
represent the computational uncertainty. The dotted line represents 
our prediction for the fractionation between olenite containing B'^' 
species only and neutral fluid. The computational eiTor is similar in 
both cases . The data points a re the values measured for dravite-fluid 
system bv lMever etaU )2008l) . 



400 - 700 °C and P = 0.2 GPa. In the experiment 
the tourmaline was represented b y dravite. In con- 



Mever et al. (2008) measured the boron isotope frac- 



tionation between tourmaline and neutral fluid at T 



trast to the former measurements of lPalmer et al.l (119921) 
the measured fractionation is very small and does not 
exceed 2.5 %c at 400 °C. Our calculated fractionation 
curve together with the experimental data are given in 
Figure |6] Our result correctly reproduces the exper- 
imental measurements within the computational accu- 
racy. The dravite-fluid fractionation is small as the 
two materials contain boron in BO3 units. We also 
predict a small fractionation between olenite carrying 
3-fold coordinated boron only and aqueous fluid, al- 
though the olenite-fluid fractionation is positive because 
of the shorter B-O bond lengths for olenite (1 .378 A vs. 
1.397 A). 



3.3.2. Boromuscovite-strongly basic fluid 

Boromuscovite s ynthesized in the experiments of 
Wunder et al.l (12005 ^ consisted of two type of polytypes, 
IM (~ 10 %) and 2M1 (~ 90 %). In order to be consis- 
tent with the experimental conditions, we derived the y6 
factors for both polytypes and computed their weighted 
average. We note that the /? factors for both polytypes of 
mica are similar with the diff'erence in (J3- 1) not larger 
than 3 %. This is consistent with similar B-O bond 
lengths (1.532 A) found in both polytypes. Boromus- 
covite contains boron in tetrahedral sites. Therefore, in 
order to investigate the impact of the B-O bond length 
on the fractionation we first compare the fractionation 
between the mineral and a strongly basic fluid contain- 
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Figure 7: The fractionation factors between boromuscovite and ba- 
sic aqueous fluid (fluid containing H4BOP. The lines represent our 
results assuming the presence of boron species in form of H4BO^ 
only (solid line) and mixture of 90% of H4BO^ and 10% of H3BO3 
(dotted line) in the fluid. The data points are the values measured by 
IWun der et al. 1 2005). The uncertainty of calculated values is indicated 
by shadowed area and is similar for both curves. 



ing boron in H4BO;j^. The result, to gether with the mea- 
surements of IWunder et al. (2005) of the fractionation 
between boromuscovite and strongly basic fluid, is sum- 
marized in Figure|7] Our calculations predict a negative 



fractionation between mica and the H4BO4 fluid. The 
agreement of our prediction with the experimental mea- 
surements is relatively good; however, the experimen- 
tal data indicate slightly stronger f ractionation. We note 
that the experimental conditions of lWunder et al.l ( l2005b 
do not assure that the measured basic fluid contained 
four-fold coordinated boron species, i.e. H4BO4, only. 
As we indicate in the Figure |7] the presence of as little 
as 10% of H3BO3 in the measured basic fluid brings the 
prediction and measurements into much better agree- 
ment. It makes sense that the (3 factor of boromuscovite 
is smaller than for aqueous H4BO4 as the average B-O 
bond length in mica is 1.532 A, while it is 1.513 A and 
therefore shorter in case of aqueous H4BO4 . 

We note that in our derivation we assumed that the 
fluid consists mostl y of ^"*^B spec i es. A s we already 



mentioned, although IWunder et alj (l2005h call the fluid 
"strongly basic" its exact composition, especially the 
amount of '^^B species is unknown. However, if for in- 
stance the '"^^B to I^^'B ratio was 1, then the predicted 
mica-basic fluid fractionation at 800 K would be 5 %c 
larger than measured. This would result in a large in- 
consistency between the computed values and the ex- 
perimental data. On the other hand, good agreement 
between the prediction and the measurements indicates 
that the strongly basic fluid was dominated by H4BO4 




1.3 1.4 
1000/T(K) 

Figure 8: The fractionation factors between boromusco vite and aque- 
ous flu id. The data points are the values measured bv lWunder et alj 
j2005t) . The solid hne represents the result obtained for ambient pres- 
sure. The dashed line represents the result for fluid containing H3BO3 
only obtained for P = 3 GPa accounting for compression and ther- 
mal expansion, with the uncertainties in calculated values indicated 
by shadowed area. The dotted lines represent the results assuming dif- 
ferent admixture of four-fold coordinated boron species (represented 
by H4BO^ with abundance indicated in the figure) to the fluid. The 
computational en'or is comparable for all the results. 



species, which is in line w ith previous studies dZeebe , 
20051: ISanchez-Valle et all 12005,) . 



3.3.3. Boromuscovite-neutral fluid 

The fractionation between boromuscovite and neu 
tral fluid involves a change in coordinat ion from B^"^J 



in boromuscovite to B^^' in neutral fluid. Wunder et al 



(I2005S measured the fractionation between the two ma- 
terials at 3 GPa, shown in Figure [8] The predicted frac- 
tionation is about 3 + 2.5 %c larger than the measured 
value. Looking for potential sources of this discrep- 
ancy, we have checked for the eff'ect of the change in lat- 
tice parameters due to combined thermal expansion and 
compression. For that purpose we applied the EOS of 
[Holland and Powell (2011) for muscovite, which gives 
4.4%, 3.8% and 3.1% decrease in volume for T - 
600 K, 800 K and 1000 K respectively and P = 3 GPa. 
As muscovites show hi ghly anisotropic compress ibil- 
ity patterns, in line with IComodi & Zanazzil d 1995b we 
applied the T and P driven change in volume assum- 
ing the 16%, 19% and 65% contribution to compres- 
sion along the a, b, ? lattice vectors. On the other hand 
(yS - 1) factors of H3BO3 and H4BO4 in aqueous solu- 
tions at P = 3 GPa increase by 3.5% and 11.8% respec- 
tively (Figures |4]and|5]l, leading to pressure-induced in- 
crease in the boromuscovite-aqueous fluid fractionation 
at the experimental pressure. The A factor, corrected 
for effects of thermal expansion and compression of 
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boromuscovite and compression of fluid, is also given 
in Figure [8] Because the high T and P eff'ects result 
in similar increases in the fi factors for both solid and 
aqueous fluid, the resulting fractionation factor between 
these two phases is close to the one derived at ambient 
conditions. Therefore, thermal expansion and compres- 
sion effects cannot explain the observed di screpancy be 



tween prediction and the measurements of iWunder et al. 
(l2005h . 

On the other hand, the comparison of our results 
with the experimental data suggests that the fractiona- 
tion between boromuscovite and fluid is the same as be- 
tween H3BO3 and H4BO4 fluids (see Figure O, which 
is at odds with the non-negligible and negative frac- 
tionation between boromuscovite and a strongly basic 
fluid. In section [3.3.21 we have shown that we are able 
to correctly reproduce the fractionation between boro- 
muscovite and strongly basic fluid, which indicates that 
our result for boromuscovite is reliable. This suggests 
that another, unaccounted effect leads to the decrease of 
the boron isotope fractionatio n between mica and neu- 
tral fluid in the experiments of IWunder et al.l ( l2005b . 

One possible solution for the discrepancy is a non- 
negligible amount of boron residing in four-fold co- 
ordinated configurations in neutral solution. This is 
in line with t he Raman spectroscopy measurements of 



Schmidt et al, (2005) , who detected a broad peak in the 
Raman spectra of neutral H3BO3 -dominated fluid and 
attributed it to B'^^ species. The integrated area of this 
peak, compared to the peak of the Raman 877 cm ' 
line of B'^*' species, indicates the presence of at least 
15 - 30 % of B'"*^ species by mole fraction. Assuming 
that there is 15 - 30 % of B'"^' species present in the fluid 
and that the p factor of these species is similar to that 
of H4BO4, the fractionation factor between boromus- 
covite and H3BO3 aqueous fluid decreases bringing the 
theory and the experiment to better agreement, which is 
illustrated in Figure |8] If this interpretation is true, it 
suggests that boron isotope fractionation could be used 
to gather information on the speciation of B in aqueous 
fluids. 

3.4. B isotope fractionation between minerals 

The boron isotope fractionation between B-bearing 
crystalline solids has received considerable atten- 



tion recently ( Wunder et al.L 2005; Meve 


r et al.L 2008 


'Klemmeetal., 201 it MarschaUl 20051; 


Hervie et al. 



.2002.) . We focus here on the investigation of boron 
isotope fractionation between mica and tourmaline as 
boron atoms in these minerals occupy sites of diff'er- 
ent coordination, which should result in a large B iso- 
tope fractionation between these two minerals. In mi- 




1.2 1.3 
1000/T(K) 

Figure 9: The fractionation factors between mica and tourma- 
line. The solid line represents our value for fractionation between B- 
muscovite and dravite. The dashed line is the experime ntal fractiona- 
tion f actor betwe en tourmaline and mica determined bv lWunder et al.l 
i2005 ) and Meye r et"ai] )20Q8|) . The experimental en'o r is 2 %o. The 
dotted line is the experimental fractionation factor of 'W under et al.l 
12005) and Meyer et al. (2008) but corrected for the presence ofBW 
species in the neutral fluid in the high P experiments of lWunder et al.l 
( 2005), as is discussed in the text. The d i amond s are the data from nat- 
ural samples taken from lKlemme et al] )201 ih and references herein. 
The uncertainties in calculated values are indicated by shadowed area. 



cas boron substitutes for silicon in the four-fold coor- 

I2OO5), while in tourmaline 



dinated site ( Wunder et al. 



(dra vite) it is incorpora ted in the three fold coordinated 
site (Meyer et al.L 120081) . Comparing the B isotope frac- 
tionation between different minerals, melts and fluids 



Wunder et al.l (120051) have shown that the fractionation 
between two materials of different B coordination is 
large, reaching 5 %c at lOOOK and much higher values 
at low er te mperatures. Klemme et al. ( 201 1 ) Marschall 
(l2005h and lHervigetal.l(l2002h measured the fraction- 
ation between coexisting phases of the two minerals 
in natural samples. The fractionation between these 
two minerals is also derived from experimenta l isotopic 



fracti onation data of B-musc ovite-fluid (Wunder et al. 
2005h and tourmaline-fluid dMever et all l2008h svs 



tems. The results of these measurements and our com- 
puted T-dependent fractionation curve are given in Fig- 
ure|9] The first striking observation is that our predicted 
fractionation factors are much larger (taking t he abso 



lute value) than the exper imental values ( Wun der et al 
2005; Mever et al. . 2008). The latter are also incon 
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sistent with the natural samples data of iKlemme et al 
( 201 1 ) and previous studies d i scusse d in that paper 
dMarschallL 120051: iHervig et all |2002l). On the other 
hand, the measurements on natural samples are consis- 
tent with our calculated values, which tends to validate 
our predictions. We notice that the most recent mea- 
surements of boron isotope signatures of tourmaline and 
white mica from the Broken Hill area in Australia by 
Klemme et alJ (1201 ll) indicate for the assumed temper- 
ature of 600 °C that the fractionation factor between the 
two phases is 10.4+2.7 %o, which is in good agreement 
with our computed value of 10.7 + 1.8 %o. The experi- 
m ental mica-tourmalin e B isotope fractionatio n factors 
of lWunder etZI (l2005h and lMever etaP (2008') are 2 %c 
and 6%o smaller at temperatures of lOOOK and 800 K 
respectively, with an experimental uncertainty of 2 %c. 
However, this discrepancy ca n be resolved by assum- 
ing that in the experiments of IWunder et alj ( 12005.) the 
fluid contained a significant admixture of B^*^ species, 
which leads to the underestimation of the experimen- 
tal boromuscovite-fluid fractionation factor by ~ 2 %c at 
1 000 K and -3.5 %c at lOOOK, as is seen in Figure [8] 
The experimental mica-tourmaline fractionation factor 
corrected for the presence of B^"^' species is also plotted 
in Figure |9l It is now more consistent with the natu- 
ral data. We note that this result independently supports 
the conclusion under lined in section [J. 3. 3 1 and result of 
Schmidt etal] (l2005h that high-f, B-beai-ing neutral flu- 
ids contain significant admixtures of B'"*^ species. 



Olenite is a mineral which can incorporate boron in 
both trigonal and tetrahedral sites as it substitutes for 
both Al and Si atoms. It is therefore interesting to check 
the fractionation of boron isotopes between the two dif- 
ferently coordinated sites in one mineral and compare it 
with the above result for mica and tourmaline. The com- 
puted fractionation between the trigonal and tetrahedral 
sites at 600 °C is 10.6 +1.9 %o, which is consistent with 
the fractionation between tourmaline and mica, indicat- 
ing that the coordination of the B atom is the driving 
factor for the fractionation of the B isotopes. Similarly, 
we computed the boron isotopes fractionation between 
trigonal and tetragonal boron sites in dravite. In order to 
create the tetragonal B site we replaced one Si atom with 
B and we added one H atom forming an additional OH 
group to compensate the charge. The computed frac- 
tionation between the sites at 600 °C is 8.9 ± 1.7 %c, 
which is also in agreement with the aforementioned re- 
sults. Next, we will show that the value of the /3 factor 
depends not only on coordination but is also strongly 
correlated with B-O bond length. 




olenite(B' '),H,BO^gas 
H,BO, fluid, dravitefs'"'') 



H^BO^' gas. H^BO^" fluid 
olenite(B'^'), dravite(B'*') 
boromuseoviie iM. 



1000 



Figure 10: /? factors for various considered materials as a function of 
temperature. The materials ai'e indicated on the right side. Their order 
reflects the value of the /3 factor at 1000 K from the largest (top) to the 
smallest (bottom). 



3.5. Fractionation between B^-^' ant/B'^^ materials 

The y6 factors computed for all the considered materi- 
als are grouped together in Figure [10] yS factors can be 
grouped into two sets, one that includes materials with 
boron in three-fold coordination and another one that in- 
cludes materials having boron in four-fold coordination. 
For olenite and dravite we also computed the /? factors 
with boron sitting on four-fold coordinated site. The p 
factor for these crystalline solids with given B'^^ /B'^'*' ra- 
tio can be derived as a weighted average of the j8 factors 
obtained for boron sitting on the two differently coor- 
dinated sites. The fractionation factor between materi- 
als of diff'erent boron coordination is ~ 8 %c on average 
at r = 1000 K. We note that it is ~ 3 %c lai-ger than 
the one deduced bv lWunder et al. ( 2005 ) from measure- 
ments performed on solids, silicate melts and fluids, 
but this can be attributed to the underestimation of the 
fractionation factors for boromuscovite-fluid system by 
Wunder et al. I (l2005h due to potential admixture of B^"^' 
species in the investigated fluid. 

Our results show also a substantial spread of fi fac- 
tors of substances containing boron of a given coor- 
dination. The spread is at least 4 %o and results from 
different B-O bond lengths. We illustrate this in Fig- 
ure [TT] by plotting together the fi factors derived for 
all considered materials at T = lOOOK as a function 
of B-O bond length. It is clearly seen that there is a 
roughly linear correlation between the p factor and B- 
O bond length, which is especially evident comparing 
the results for crystalline solids. For instance, out of the 
considered B'^^^-bearing minerals boromuscovite has the 
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Figure 11: The /j factor at T = lOOOK for various considered ma- 
terials as a function of B-O bond length. Filled circles represent the 
values obtained for boron in trigonal sites and open circles represent 
the values obtained for tetragonal sites. 



longest B-O bond length of 1.525 A (IM) and 1.516 A 
(2M1), followed by dravite 1 .5 14 A and olenite with B- 
O bond length of 1 .502 A. This tracks the differences in 
the /3 factors derived for these materials. In addition the 
materials having B'^^'' species only exhibit shorter bond 
lengths of ~ 1 .37 A and higher /3 factors, while the ma- 
terials containing B'^^ species having bond lengths of 
about ~ 1.52 A show much smaller /3 factors. There- 
fore, the tighter bonding of B'^' species likely explains 
why the heavy B isotope prefers the less coordinated 
phases. This clearly shows that the change in the B-O 
bond length during an isotope exchange is the leading 
factor driving the production of the boron equilibrium 
isotope signatures at high T. 



4. Conclusions 



Bigeleisen & Maven (1 1947b "single atom approxima- 
tion " method we demonstrated that using the pseudofre- 
quencies derived from the force constants acting on the 
fractionating element together with the full formula for 
computation of the reduced partition function ratios re- 
sults in significant improvement in the accuracy of the 
computed fractionation factors, which is essential when 
lower temperature materials and high vibrational fre- 
quency complexes are considered. 

In order to understand the fractionation between B- 
bearing crystalline solids and aqueous fluids we per- 
formed a set of calculations of p factors for dravite, 
olenite, boromuscovite and aqueous solutions of H3BO3 
and H4BO4. In agreement with the experimental find- 
ings we show that the fractionation strongly correlates 
with coordination through the change in the B-O bond 
length. The lower trigonal coordination BO3 arrange- 
ment results in higher "B/'^B (by ~ 8 %c at T = 
lOOOK) than the tetrahedrally coordinated boron com- 
plexes, which exhibit ~ 0.15 A longer B-O bonds. The 
computed fractionation between minerals and fluids of 
the same coordination are in good agreement with ex- 
periments. However, we predict larger isotope fraction- 
ation between boromuscovite and H3BO3 fluid (by at 
least a few %o) than was measured in situ at high P by 



Wunder et al.l (,2005) and lMever et al.1 (2008), but that is 
consistent with measurements on natural samples. We 
note that the presence of B'^^ in high-P fluid could rec- 
oncile the in situ experimental results with our predic- 
tion and other r neasurements. This is expected from the 
experiments of 'Schmidt et alj (2005), but requires fur- 
ther experimental confirmation. If true, this would open 
the possibility for using the isotope fractionation tech- 
niques as a tool to measure the speciation of boron in 
fluids and crystalline solids. We have also demonstrated 
that with our computational approach we are able to cor- 
rectly predict the pressure-induced isotope fractionation 
for compressed aqueous fluids, which indicates the abil- 
ity of ab initio methods to predict the isotopic signatures 
of highly compressed materials, even those that are dif- 
ficult to investigate experimentally. 



In this work we have presented a detailed anal- 
ysis of boron isotope fractionation between boron- 
bearing crystalline solids and aqueous fluids at high 
T and P conditions. In order to perform our in- 
vestigation we have applied and extended a computa- 
tionally efficient approach for the computation of iso- 
tope fractionation factors for complex minerals and 
flu ids at high temp e rature s and pressures presented 
by Kowalski & JahnI ( 201 lb . As an extension to the 



Our study confirms that ab initio computer simula- 
tions are a useful tool not only for prediction but also 
understanding the equilibrium stable isotope fractiona- 
tion processes between various phases, including aque- 
ous solutions, at high pressures and temperatures. They 
can nicely complement experimental efforts, provide 
unique insight into the isotope fractionation process on 
the atomic scale and deliver data for conditions that are 
inaccessible by the current experimental techniques. 
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